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ABSTRACT 

We propose two detection techniques that take advantage of a small sky area approximation and are based on modifi- 
cations of the internal linear combination (ILC) method, an approach widely used in Cosmology for the separation of 
the various components that contribute to the microwave background. The main advantage of the proposed approach, 
especially in handling multi-frequency maps of the same region, is that it does not require the a priori knowledge of 
the spatial power-spectrum of either the CMB and/or the Galactic foreground. Hence, it is more robust, easier and 
more intuitive to use. The performance of the proposed algorithms is tested with numerical experiments that mimic the 
physical scenario expected for high Galactic latitude observations with the Atacama Large Millimeter/submillimeter 
Array (ALMA). 

Key words. Methods: data analysis - Methods: statistical - Cosmology: cosmic microwave background 



1. Introduction 

The detection of (extragalactic) point-sources in experi- 
mental microwave maps is a critical step in the analysis 
of the Cosmic Microwave Background (CMB) maps. Beside 
the specific interest related to the construction of dedicated 
catalogues, if not properly removed, these sources can have 
adverse effects on the estimation of the power-spectrum 
and/or the test of Gaussianity of the CMB component. For 
this reason, this subject has been extensively considered in 
literature. Various detection techniques have been proposed 
but no general consensus about their real performances and 
properties has been reached. The reason is due to the dif- 
ferent approaches adopted for dealing with the CMB and 
foregrounds (both extragalactic and Galactic) that make 
difficult the development of a general theory. Often such 
contributions are considered as a "noise" to be added to 
the instrumental noise. The various methods differ in the 
statistical characteristics attributed to such noise as well as 
in the way to deal with it. 

Many efforts have been dedicated to the case of multi- 
ple frequencies maps of the sam e sky-area and many algo- 



rithms have be en proposed (see Herranz and Sanz 
2012 



Herranz et al. 



j2QQ8 

and re ferences therein). Apar t from 



a recent Bayesian approach ( [Carvalho et al. 



|2QQ9[ ), most 

of them belong to two broad classes of techniques. The 



first class makes use on the Ney man- Pear son (NP) crite- 
rion that consists of the maximisation of the probability of 
detection Pd under the constraint that the probability of 
false alarm PpA (i-e., the probability of a false detection) 
does not exceed a fixed value a ( |Kay ||1998[ ). The result- 
ing algorithms are an extension of the classic matched filter 
(ME). The second class is based on the constrained max- 
imisation of the signal-to-noise ratid^ (SNR) of the source 
intensity with respect to the underlying background. The 
constraints are chosen in such a way to improve detection. 
This class provides algorithms of ''internal linear combina- 
tion^^ (ILC) type. In Cosmology ILC is essentially used for 
the separation of the various components contributing to 
the microwave sky emission (Eriksen et al'~||2QQ4 Hinshaw 



et al. ||2QQ7[ |Vio and Andreani. ||2QQ8p . The main limita- 



tion of both classes is the requirement that the background 
is made of realisations of stationary stochastic processes 
with spatial spectra that, in addition, are supposed to be 
known. For example, on small patches of sky, the Galactic 
contribution is modelled with stationary, in the case of NP 
algorithms Gaussian, processes with steep spectra (e.g. 1// 
noise). These are rather rough assumptions. For this reason, 
in the context of futu re high spatial resolution observations 
Ramos et al. (2011) propose a method for high Galactic 
observations where the maps, almost free of the Galactic 
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contamination, are linearly combined in such a way that 
the resulting one is free of the CMB emission. Although 
such method has been successfully applied to the WMAP 
maps, it suffers from the drawback that it is unsuited for 
the detection of point-sources with spectrum similar to that 
of the CMB or the Sunyaev-Zeldovich effects. In this paper, 
we follow a different approach, based on two algorithms to 
be used in conjunction, that does not require the a priori 
knowledge of the spatial power spectra of the diffuse emis- 
sions due to either CMB or to the Galaxy and at the same 
time it permits to deal with the above mentioned prob- 
lem. The basic consideration is that the detection of point 
sources is typically done on very small areas of sky where 
both the CMB and the Galactic components can be very 
well approximated with low-degree two-dimensional poly- 
nomials. This is a more reasonable assumption than the 
ones mentioned above, i.e. emissions resulted from realisa- 
tion of stochastic spatial processes. The performances of 
the algorithms is tested via numerical experiments based 
on simulated maps of high Galactic latitude that might be 
the area of interest of CMB high spatial resolution obser- 
vations. 



2. Formalization and solution of the problem 

When looking for a point source in an area of sky, data 
can be modelled as two-dimensional discrete maps 
each of them containing Np = Np^x Np^ pixels, correspond- 
ing to Nf different observing frequencies (channels), with 
the form 

Xi = Si^Ci^Zi^gi^Ari. (1) 

Here, Si is the contribution of the point-sources at the ith 
frequency, C^, and Gi are the backgrounds due to CMB, 
extragalactic and Galactic emissions, respectively, and Afi 
is the instrumental noise. In this model, the contribution of 
the point-sources is assumed in the form 

Si = aiT, (2) 

with "a^" the intensity of the source at the zth channel. 
According to Eq. ([2| , and without loss of generality, all the 
sources are assumed to have the same profile T indepen- 
dently of the observing frequency. In practical applications, 
this is not true. However, it is possible to meet this condi- 
tion by convolving the images with an appropriate kernel 
with no remarkable consequences (see below). 

For computational reasons, that soon will become evi- 
dent, it is useful to convert the two-dimensional model ([T]) 
into the one-dimensional form 

Xi = Si^ Ci^ Zi^ g^^ Ui. (3) 

Here, Xi = VEC[Ari], with VEC[i^] the operator that 
transforms a matrix H. into a vector by stacking its columns 
one underneath the other. Something similar holds for the 
other quantities. 



2.1. Detection with ILC background removal 

The main issue in detection problem when more maps of the 
same area are available at different observing frequencies, 
is how to handle them. One classical solution consists to 
linearly compose the maps by means a set of weights w = 



with a single map given by 



^ in such a way that it is possible to work 



X 



(4) 



where X = [xi, X2, . . . , xat^] is di Np x Nf a matrix. The 
question is how to fix such weights. In accomplishing such 
a task, it is necessary to take into account that there is 
some a priori information about the various components in 
Eq. ([T]). In particular: 

— For each observing frequency z, the spectra of Ci and 
3i are known with good accuracy. This is the case for 
instance of the CMB and the Sunyaev-Zeldovich (SZ) 
effects, both thermal and kinetic; 

— Ci and Gi have a diffuse spatial distribution with a 
spatial scale much greater than that of the point sources; 

— Noises Afi can be reasonably assumed as given by the re- 
alisation of independent Gaussian white-noise processes 
with standard deviation o-/^.. 

We are interested in exploring the situation in which 
in addition to the CMB the extragalactic component con- 
sists of secondary anisotropics of the CMB. Here we con- 
sider only the SZ effects which are the strongest ones in 
gala xy clusters, groups of galaxies and in protoclusters 
(Birkinshaw [1999 ) and whose spectral shape is well known. 
The first point implies that the contribution 6 to of the 
CMB and of the SZ components can be obtained from 



b = M[c,z], 



(5) 



where c and z are templates for CMB and SZ components 
(i.e. maps that do not depend on frequency) and M is a 
Np X Ne mixing matrix. In the present context, TVg = 2 
since the kinetic SZ emission has the same spectrum of 
CMB. For this reason, from now on, with Ci we will indicate 
the CMB plus the kinetic SZ emission. The second point 
implies that within a small area centred at a point source 
the CMB and the Galactic emissions vary very little. This 
suggests that, for any sub-map ^i{j^ k) with —Nj < j < Nj 
and -Nk < k < Nj, {Nj <C Np^ and Nk <C A/pJ, these 
emissions can be safely approximated by a low-degree, two- 
dimensional polynomial of degree m 



'Pm{j,k)=Y,Mfkl; q + r<l, 



(6) 



^=0 



where {ai} are real coefficients whereas q and r are integer 
numbers permuted accordingly. 

If for detection one adopts the criterion of the SNR max- 
imisation, all these considerations suggest a model where 
the weights, the point source intensity and the parameters 
of the approximating two-dimensional polynomial in x are 
optimised simultaneously, i.e. Q 



R{w^ c. A) 



arg mm 

w.c.X 



(Xw — Lc) 



^X^{[a,Mfw-ei] 
(7) 

"argminF(x)" and 



^ We recall that the functions 
"arg maxF(x)" provide the values of x for which the function 

X 

F{x) has the smallest and greatest value, respectively. 
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Here, "||.||" indicates Euclidean norm, a = 
[ai, a2, . . . , ttAT^]-^ is an array containing the emission 
spectrum of the source to detect, c = [a,a^]^ is an 
array with size Nc = [{m + l)(m + 2)/2] + 1 with 
"a" the amphtude of the source in the hnearly com- 
posed sub-map of AT and a the coefficients of the 
two-dimensional polynomial. L is a Np x Nc matrix 
with the form L = [/, P] where / = VEC[:r] and P 
the Np X {Nc — 1) matrix that results from the least- 
square fit of a two-dimensional polynomial. For example, 
when m = 1, P = [^1,^27!] whereas, when m = 2, 
P = [Si 061,62 062,61 0^2,^1,^2,1], where "0" repre- 
sents the element-wise matrix multiplication (Hadamard 
product), 1 is a vector of ones and = VEC[Ai], 
^2 = VEC[A2] where Ai is a matrix with 2Nj + 1 

identical columns [— A^fe, —J^k + 1, • • • , 0, . . . , A^/c — 1, J^k]^ 
whereas A2 is a matrix with 2Nk + 1 identical rows 
[-Nj, -Nj + 1, . . . , 0, . . . , AT,- - 1, A/"^-]. Finally A is a A^e + 1 
array of Lagrange multipliers whereas ei is a A'g + 1 array 
of zeros except for the first element that is "1". With 
this model the quantity SNR = {a^w)'^/\\{Xw - Lc)\\'^ 
is maximized under the constraint that a^w = 1 (i.e., 
preservation of the source amplitude) and M^w = [0, 0] 
(i.e., the CMB and SZ components are zeroed). 

The basic idea behind this method, that we call modi- 
fied multiple ILC (MMILC), is that, if in the centre of the 
selected sub- map there is a point-source, then the value of 
"a" should exceed a threshold due to noise. After some al- 
gebra, one obtains that the solution of problem ^ is given 
by the system of equations 

+2Cxx -2CxL /w\ /0\ 

-2CxL ^2Cll c = , (8) 

M J \ Xj \eiJ 

where Cxx = X^X, Cxl = X^ L and Cll = L^L. The 
explicit solution for c and A is not difficult to obtain 
but it produces rather complicated expressions. Hence, the 
numerical solution is more advantageous. One interesting 
characteristic of solution ([8| is that it does not require the 
knowledge of the noise level of each map, a quantity that 
often can be only roughly estimated. 

When searching for point sources in a given set of maps, 
the procedure consists in fixing the size {2Nj^l) x (2A^/e + l) 
of a window that is made to slide, pixel by pixel, across the 
area of interest. At the end of this procedure a single map 
is obtained containing the estimated values of "a" for each 
pixel. Now, the question is to fix the detection threshold 
below which a given value of "a" is supposed to be due 
only to noise. In this respect, the direct use of solution (|8| 
is difficult. For this reason, two different procedures are 
suggested: 

1. a = if a < /cdL, where /c is a constant fac- 
tor (typically k = 4,5), (Jl = \\alw\\^J {L^ L)'^^, 
o-n = [crni,crn2, . . . ,crn^J^ and {L^ L)^^^ is the first 

entry of matrix {L^L)~^. This operation corresponds 
to estimate the standard deviation da of "a" for a fixed 
It;. Such an approach has the advantage that the matrix 
{L^ L)~^ can be computed only once since it is the 
same for all the sub-maps. But, it has the disadvantage 
that the standard deviations of the noises {n^} are to 



be known in advance; 

2. a = if a < /camap, where again A: is a constant factor 
and (Jmap is the standard deviation of the entries in the 
final map. This is an unsophisticated approach, however 
it has the advantage that does not require the standard 
deviation of the noise in each sub-map, a quantity usu- 
ally only roughly known. 

Perhaps an advisable procedure consists in using both 
methods and to check for differences. 

2.2. Detection without ILC background removal 

The MMILC detection procedure is potentially quite effec- 
tive, however it suffers of two main drawbacks: 

1. In order to remove the CMB and SZ components, one 
or more of the weights in w have to be negative. As 
a consequence, since in the final map a = a^w and 
<^map = ll^n'^ll^ givcu by the sum of positive as 
well negative values whereas (Jmap is given by the sum 
of positive values only. In other words, the background 
subtraction reduces the SNR with respect to a simple 
sum of the maps. The situation worsens when the 
emission of a point source has a spectrum similar to 
that of CMB or of SZ since "a" will tend to zero; 

2. If a is an array such that Ma = 0, i.e. a belongs to the 
nullspace of M (i.e., it is given by the linear combina- 
tion of the column of M) then the system ^ does not 
have any solution because the constraints a^w = 1 and 
M^w = become incompatible. 

For this reason, in order to detect point sources with a be- 
longing to the nullspace of Af , the above procedure has to 
be adapted to work without the ILC removal of the CMB 
and SZ components. More specifically, problem ^ is con- 
verted into 

R{w,c, A) = argmin [\\{Xw - Lc)\\'^ + X{a^w - 1)] , 

w,c,X 

(9) 

with solution given by 

/ +2Cxx -2CxL a\ /w\ /0\ 
-2CxL +2Cll c = . (10) 

V 0^ ; V A / V 1 / 

With this method, that we call modified ILC (MILC), the 
CMB is not removed through the use of the mixing matrix 
M. However the fact that this is a component with diffuse 
spatial distribution makes us hope that it can be removed 
through the polynomial approximation of the background. 
As a consequence, in the final map the only contribution 
beyond that of the point sources is the SZ (both thermal 
and kinetic) . Unfortunately this is an unavoidable problem. 
Without further information it is impossible to separate 
an SZ emission with point like spatial distribution from a 
genuine point source. 

3. Practical uses 

In this section we discuss some practical problems and how 
they can be addressed. The first is related to the degree 
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m of the polynomial used to approximate the background. 
We need to take into account that, when two polynomials 
of degree mi and m2 with mi > m2 are summed together, 
a polynomial with degree m = mi is obtained. Hence, m 
is fixed by the diffuse component that requires the largest 
degree of the approximating polynomial. If the sky patches 
are small, it can be reasonably expected that a first degree 
polynomial is a good choice. The second question is related 
to the sizes Nj and Nj^ of the sub-map where to test for the 
presence of a point source. Two competing requirements 
raise: on the one hand Nj and must be as large as pos- 
sible to reduce errors in the estimation of the polynomial 
parameters, on the other hand, a small size implies that 
the polynomial approximation for the background is pre- 
ciser and a lower probability is expected with respect to 
the area where some other sources are present. For illus- 
trative purposes. Fig. [l] shows the standard deviation aa of 
the estimated intensity a as provided by MILC in the case 
of a point source with a Gaussian profile and a dispersion 
cTpsf equal to 3 pixels. A single map is considered where 
the background is given by a two-dimensional one degree 
polynomial, a Gaussian and white instrumental noise with 
standard deviation (Jn, and A^^^, Np^ are increased. The 
true value of a is 1 in unit of The decrease of a a is evi- 
dent. Figure [2] shows the relationship between Pd and PpA 
for different values of the ratio a /an- From these figures it 
is clearly shown that Nj, Nj^ lying in the range 3crpsf-5(Tpsf 
is a reasonable compromise. 

Another issue is related to the fact that in practical ap- 
plications the width of the PFSs changes with observing 
frequency. Widespread practice is to convolve maps with a 
suited kernel function in order to get a common spatial res- 
olution at all frequencies. This operation has the beneficial 
effect to reduce the standard deviation of the instrumental 
noise but at the same time it introduces a spurious spatial 
correlation in it. Actually, even if neglected, this latter is 
not critical since MILC and MMILC are linear techniques 
and the only consequence is a slight reduction of their effi- 
ciency. In other words, given the above mentioned reduction 
of the standard deviation of the noise this spurious corre- 
lation is of secondary importance. This is especially true if 
one takes into account that there are other and more im- 
portant approximations that make the analysis of data less 
rigorous (e.g., often the level of instrumental noise is only 
roughly known). 

As a final comment, both MMILC and MILC work opti- 
mally only for a specific emission spectrum a, a limitation 
common to many other detection techniques. But this is 
not critical. It is sufficient to apply the detection algorithm 
each time with a different value of a. This is made possible 
by the fact that MMILC and MILC are quite fast algo- 
rithms since they require the numerical solution of systems 
containing no more than a couple of tens of linear equa- 
tions. 



4. Numerical experiments 

In order to support the arguments presented above, here we 
present some numerical experiments with simulated maps 
at high Galactic latitude (where the Galactic contamination 
is negligible) that is the region of interest for future CMB 
experiments. 

We produced small sky patches of 0.86 deg^ at 3'^ angu- 
lar resolution with several components, namely, the CMB 



and the Sunyaev-Zel'dovich effects (SZ), both kinetic and 
thermal. To produce these maps we used Hydrodynamic/N- 
body simulations with cosmological parameters consistent 
with WMAP parameters for a flat Universe and standard 
A CDM model, with an equation of state for the dark en- 
ergy component of w = —1. The adopted present time den- 
sity parameters expressed in terms of the critical density 
are (^^cdm, ^a, ^5) = (0.256,0.7,0.044), a dimensionless 
Hubble constant of h = 0.71 and a mean CMB temperature 
of T=2.725 K. It is assumed adiabatic initial conditions, a 
spectral index of = 1 and full reionisation at redshift 7. 
For the present epoch we considered a normalisation power 
spectrum of ag = 0.9 and a shape parameter of F = 0.17. 
The CMB compon ent is produced with the CAMB code 



(Lewis et al. 



trum. The fn. 



i 



2000|) to obtain the linear CMB power spec- 
was 



sky CMB temperature anisotropy map 
generated with the HEALPix software (Gorski et al. 12005) 



with Nside = 8192. From this map it was extracted a small 
sky region with an area of about 0.86 deg^ around the equa- 
tor, projected in a squared map. Details about th e simula- 
tions of the SZ effect components can be found in (da Silva 



|et al. ||2001| [Ramos et al. ||2012p . The frequencies chosen 
were 90, 150, 250, 330, 440, 675 and 950 GHz which cor- 
respond to the receiver bands of ALMA. All components 
were co-added resulting in a final map, A/cmb+sz/^, with 
pixel size of 3 arcsec. We use the central part of the maps 
(300 X 300 pixels) and convolve them for each frequency 
with a Gaussian PSF with dispersion of 3 pixels. To each 
map a white-noise process with standard deviations (7^ set 
to 0.12 time the standard deviation of the values of map 
itself has been also added. Finally 20 point sources ran- 
domly distributed have been included with = 1.7 a m- 
In this way, maps with the same SNR are obtained. The 
values of cr^^ and have been arbitrarily chosen in such 
a way to test algorithms under very bad operational con- 
ditions but at the same time to obtain stable results (i.e. 
with different realizations of the noise process almost all the 
sources are correctly detected with no false detections). The 
simulated experimental scenario corresponds to an adverse 
situation of rather low SNR and, since a^. increases with 
frequency, with a spectrum a that mimics that of the CMB 
plus SZ background (le. a is close to the nullspace of M, 
or Ma ~ 0). Figures ISp] show the noise free and the noisy 
version of such maps, respectively. From these figures it is 
evident that most of the point sources are not even visible 
and that are by far exceeded by the SZ point-like emis- 
sion. Figure |5] shows the results obtained with MILC and 
MMILC using a detection threshold 5crL and a background 
approximated by a two-dimensional first degree polynomial. 
As expected, the MMILC does not work. On the other 
side, MILC has effectively removed the CMB components 
and correctly detected all the point-sources in the map. 
However, many of the SZ point-like emission is also present. 
The situation greatly improves when the value of amplitude 
of the point sources at 90 GHz is set to zero. In this way 
Ma 7^ 0, i.e. the degeneracy of the similarity of CMB and 
SZ kinetic spectra is broken. Similar results are obtained 
setting the flux to zero at another frequency. As visible in 
Fig. |6j while MILC provides almost the same results as be- 
fore, MMILC is able to correctly detect all the point sources 
as well to completely remove the CMB and SZ contamina- 
tion. A point to stress is that, in the present high Galactic 
latitude experiment the diffuse Galactic component is set 
to zero. Hence, a simple ILC algorithm could be used (i.e. 
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without the subtraction of the polynomial approximation 
of the background). However, the use of MMILC permits 
to test this algorithm for higher level of noise. 

5. Summary and conclusions 

In this paper two different modifications of the internal lin- 
ear combination algorithm, the modified ILC (MILC) and 
the modified multiple ILC (MMILC), have been presented 
to detect point sources in multifrequency CMB maps. The 
remarkable property of the proposed algorithms is that they 
do not require the a priori knowledge of the statistical 
characteristics of the diffuse sky background since this is 
assumed to be locally approximate to a low-degree two di- 
mensional polynomial. In fact, the methods available in lit- 
erature have the drawback to assume this component as the 
realisation of stationary, often Gaussian, two-dimensional 
stochastic processes, with known spatial power-spectrum, 
which is not a realistic approximation. The two proposed 
algorithms used in conjunction are effective in the detection 
of point sources independently of the spectral characteris- 
tics of their emission. Their potential performance has been 
tested via numerical experiments. 
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Fig. 1. Standard deviation aa of the estimated intensity a as provided by MILC in the case of a point source, with 
Gaussian profile with dispersion (jpgf equal to 3 pixels, as a function of the sizes Nj = Nk of the searching sub-map. Here, 
a single map is considered with a background given by a two-dimensional one degree polynomial, instrumental noise is 
Gaussian and white with standard deviation cr„. The true value of "a" is 1 in unit of (7„. 




Fig. 2. Relationship between the probability of detection^ Pd, vs. the probability of false alarm, PpA for the case shown in 
Fig. fl] but for different values of the ratio a/an- Note the different scale used for the abscissa in the bottom-right panel. 
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Fig. 3. Noise free maps simulating a high Galactic declination area of sky at the ALMA observing frequencies. 20 
randomly distributed point sources with the same intensity have been added. The point source intensity has been set to 
1.7 times the standard deviation of the noise (see next figure). In this way their spectrum mimics that of CMB + SZ 
(see text). The PSFs are assumed to be Gaussian with a standard deviation of 3 pixels. The two bottom-right panels 
show the simulated point sources and their position on the 950 GHz map, respectively. 
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Fig. 4. Noisy version of the maps in Fig.js] Noise is Gaussian- white with standard deviation set to 0.12 time the standard 
deviation of the values in the corresponding noise free maps. Ah of the point sources have the same intensity set to 1.7 
times the standard deviation of the noise. The two bottom-right panels show the simulated point sources and their 
position on the 950 GHz map, respectively. 
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Detection map — MMILC Detection map — IVIILC 




Fig. 5. Results provided by MILC and MMILC when applied to the maps in Fig. [4j The detection threshold has been 
set to SdL (see text) and the background has been approximated with a two-dimensional polynomial of degree one. The 
bottom left panel shows clearly that MMILC is not able to retrieve point sources in this case. This happens because 
the point surges spectra has a frequency-dependence similar to that of the CMB and SZ and therefore the subtraction 
process gets rid of all of them. The bottom right panel shows that MILC, on the contrary, retrieve all sources and the 
SZ point-like emissions, because it subtracts the underlying diffuse component with the polynomial approximation. 
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Fig. 6. As in Fig. |5| with the only difference that the intensity of the sources at 90 GHz is set to zero. With this constrain 
both MMILC and MILC are able to detect all the point sources and get rid of the SZ point-like emissions. 



